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Abstract 

Local density approximation (LDA) to the density functional theory (DFT) has continuous 
derivative of total energy as a number of electrons function and continuous exchange-correlation 
potential, while in exact DFT both should be discontinuous as number of electrons goes through 
an integer value. We propose orbital densities functional (ODF) (with orbitals defined as Wannier 
functions) that by construction obeys this discontinuity condition. By its variation one-electron 
equations are obtained with potential in the form of projection operator. The operator increases 
a separation between occupied and empty bands thus curing LDA deficiency of energy gap value 
systematic underestimation. Orbital densities functional minimization gives ground state orbital 
and total electron densities. The ODF expression for the energy of orbital densities fluctuations 
around the ground state values defines ODF fluctuation Hamiltonian that allows to treat correlation 
effects. Dynamical mean-field theory (DMFT) was used to solve this Hamiltonian with quantum 
Monte Carlo (QMC) method for effective impurity problem. We have applied ODF method to the 
problem of metal-insulator transition in lanthanum trihydride LaH3_a;. In LDA calculations ground 
state of this material is metallic for all values of hydrogen nonstoichiometry x while experimentally 
the system is insulating for x < 0.3. ODF method gave paramagnetic insulator solution for LaHa 
and LaH2,75 but metallic state for LaH2,5. 

PACS numbers: 78.70.Dm, 71.25.Tn 
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I. INTRODUCTION 



Numerical electronic structure calculations are now well established branch of solid state 
physics. While for the finite systems, such as atoms and molecules, more sophisticated 
and rigorous calculation methods exist, for extended systems studied in condensed matter 
physics the only widely used practical tool is up to now density functional theory (DFT) 
in the local density approximation (LDA)i^ It has so great predictive power, that not only 
charge and spin density, one-electron and total energies obtained in LDA are in general 
in a very good agreement with experimental data, but it was possible to develop ab initio 
molecular dynamic methods, based on LDA. Such methods achieved the level of numerical 
experiment, because even so complicated effects, such as reconstruction of the crystal surface 
can be correctly described by them.^^ 

However, there are materials where LDA results do not agree well with experimental data. 
For band insulators and semiconductors LDA gives systematically underestimated values of 
the energy gap.— For Mott insulators, for example transition-metal oxides, LDA results could 
be qualitatively wrong, giving metallic state while experimentally those systems are wide 
gap insulators.— 

There were many attempts to cure this deficiency of LDA. Among the most widely used 
one can mention GW,— SIC,— and LDA+U^ methods. While those approaches have their 
advantages there is still no universally accepted calculation scheme which would be as sim- 
ple and practical as standard LDA and a search for better methods continues in scientific 
community. 

The basic problems of LDA can be traced to the fact that the exchange-correlation energy 
functional is defined in a local approximation. As a result its variational derivative, exchange- 
correlation potential (that is a function of the density value in particular point r instead 
of being a general functional of density), is continuous as a number of electrons function. 
More general approximations using in addition to the electron density its gradient like GGAi^i 
method have continuous potential as well. However, Perdew et al^ had investigated the 
properties of the exact density functional (EDF) and shown that its potential must jump 
discontinuously when number of electrons goes through an integer value. The proper 
function of total energy E versus number of electrons should have a curve as series of 
straight-line segments with derivative discontinuities at integral values of while in LDA 
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this curve has continuous derivatives. Any attempts to improve LDA as an approximation 
to exact density functional theory should be done in such a way that a new functional would 
obey this discontinuity requirement. 

In the present work we define a functional that by construction has the potential discon- 
tinuities required for exact density functional. For that we introduce a concept of orbital 
densities corresponding to one-electron orbitals. The orbital densities functional (ODF) 
depends on a set of orbital densities instead of the total electron density only. When the 
number of electrons goes through an integer value the variational derivative of the ODF func- 
tional jumps discontinuously and the corresponding function of total energy Eqdf versus 
number of electrons has a curve as series of straight-line segments. 

Varying this functional one-electron equations were obtained with potential in the form of 
projection operator. For occupied states the effect of this operator on the energy is negative 
and positive for empty states. As a the result it increases separation between valence and 
conduction bands comparing with LDA values thus curing LDA deficiency of energy gap 
value systematic underestimation. 

The minimization of ODF functional gives a set of orbital densities and hence the total 
electron density corresponding to ground state of the system. The same functional can be 
used to calculate energies of orbital densities fluctuations around the ground state values. 
From this we have derived fluctuation Hamiltonian defined via orbital density fluctuation 
operators. This Hamiltonian allows to treat correlation effects and hence obtain better 
description of ground state properties and spectral function for the excitations. In the present 
work we have used for ODF Hamiltonian dynamical mean-filed theory (DMFTy^*i^^ with 
quantum Monte Carlo (QMC) method to solve the effective impurity problem. 

The optimal choice for one-electron orbitals needed to define orbital densities can be 
determined by condition of fluctuation energy minimum. We have shown that the less 
extended in space these orbitals are, the lower the energy of orbital densities fluctuations 
around ground state is. Therefore in our ODF method we have used maximally localized 
Wannier functions for orbital densities definition.i^ 

Recently we have developed "generalized transition state" (GTS) methocU^ to improve 
agreement of calculated and experimental spectral properties comparing with LDA. We have 
found that ODF projection operator is identical to the one in equations of GTS method. 
The basis of GTS method was an idea that one-electron energies corresponding to Wannier 
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functions should have a meaning of the removal (addition) energies for electrons from (to) 
the corresponding states. That was realized by using "transition state" schemo^^ gener- 
alized to Wannier functions. It is remarkable that so different approaches as "transition 
state" correction to excitation energies and requirement of exchange-correlation potential 
discontinuity can lead to the same equations. 

We have applied ODF method to the problem of metal-insulator transition in LaHs-a;. 
LDA has severe difficulties for this material because it gives metallic solution for all values 
of hydrogen deficiency x while experimentally the system is insulating for x < 0.3. We have 
found that ODF projection operator potential is enough to open a gap for stoichiometric 
LaHa but in order to reproduce paramagnetic insulator for LaH2.75 correlation effects should 
be taken into account via DMFT-QMC solution of ODF fluctuation Hamiltonian. 

The paper is structured as follows. In Sec. |n] the derivations of ODF functional and 
fluctuation Hamiltonian are presented. In Sec. IIIII and IIVI we describe the construction of 
Wannier functions. Sec. |3 gives a short account of DMFT method and Sec. IVII describes 
the self-consistency issues. In Sec. IVIII results of ODF calculations for LaH3_a; are presented. 
Sec. IVIIII concludes the paper. 

II. ORBITAL DENSITIES FUNCTIONAL AND FLUCTUATION HAMILTO- 
NIAN 

In the local density approximation (LDA)^ to the density functional theoryi(see Ap- 
pendix [X]) exchange-correlation energy is calculated via: 



where £a;c(p(i")) is exchange- correlation energy density for homogeneous electron gas with 
density equal to p(r). 

Equation ((H) deflnes an exchange-correlation potential V^c(p(r)) = 6{Exc[p])/Sp{r) which 
is a continuous function of the number of electrons N. However, Perdew et al^ proved that 
for the exact density functional exchange-correlation potential must jump disco ntinuously 
when number of electrons goes through an integer value. 

The Hohenberg-Kohn theorent^i was extended in Ref. Q to fractional electron number. 
It was shown that for the electron density p(r) which integrates to = M + u; where M is 




(1) 
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an integer and < c<j < 1 the exact density functional Eedf is equal to: 

Eedf{M + uj) = {1- uj)Eedf{M) + ujEedf{M + 1). (2) 

This means that in general the curve of Eedf versus is a series of straight-line segments 
with derivative discontinuities at integer values of A^. 

It was provenii that the chemical potential /i = OEedf /dN is discontinues when the 
number of electrons goes through integer value: 

r_/ (M-l<iV<M) 

= < (3) 

y-A {M <N <M + 1) 

where / = Eedf{M - 1) - Eedf{M) and A = Eedf{M) - Eedf{M + 1) are removal and 
addition energies respectively. The functional derivative 6Eedf/Sp{t^) is also discontinuous: 
two limits for approaching M from above and below will differ by constant I — A. 

Electron density [Eq. ()A1|) ] can be expressed as a sum of the "orbital densities" Pi{r) 
defined in the following way (ipiir) is i-th orbital wave function): 

p(r) = J]p,(r), (4) 

i 

Pi{r) = ni|^i(r)|^. 

The orbital densities can be varied from zero to the maximum values p™'^^(r) = |?/;j(r)|^. 

The condition of linear dependence of the exact density functional on the fractional num- 
ber of electrons [Eq. can be expressed via orbital density Pj{r) corresponding to the 
partially occupied orbital j. This orbital is the lowest unoccupied one for number of elec- 
trons N = M and the highest occupied for N = M + 1. The variation of the total electron 
density p(r) will be defined only by the variation of Pj{r) for number of electrons changing 
from M to M + 1 

Eedf[p] = Eedf[po + Pj] = Eedf + / drpj{r)Aj{r), (5) 

PjW=o J 

where po(r) is the density for number of electrons A^ = M. Here Aj(r) satisfies to the 
equation: 

/ drp7'^^(r)A,(r) = Eedf - Eedf • (6) 
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Orbital density pj (r) enters Eq. (jSJ in a linear form. However LDA functional [Eqs. ()A2|1 
and (Q)] does not show such a linear dependence on density variations. Equation analogous 
to Eq. (0) for LDA has a general form (keeping only first and second variational derivatives 
in expansion series): 

, 5 Eld A 



Elda[p\ = El£)a[po + Pj] ~ Elba 



p,(r)=0 



Spiv) 

\ J drpjir) j dv'pjir') ^ ^^^^ 



Pj{r)=0 
p,{r)=0' 



(7) 



(5p(r)5p(r') 

Using Eqs. ()A24|) and ()A25|) this expansion can be rewritten via constrain potential $(r) 
and effective interaction strength function W{v^v'): 



Elda[p\ = EldaIpo + Pj] ~ Elda + / rfrpj(r)$(r) 

Pi(r)=0 J 

+ ^ J dvp,{r) j dr'p,{v')W{ry) 
If one will add to Eq. ((7j) a correction term: 

Ecorr[pA^-\ j drp,{v) j dv'{p,{r')-p^'^%r'))W{vy, 



Pj(r)=0 
p,(r)=0' 



(8) 



p,(r)=0 



(9) 



then the linear dependence on the density variations (as it is required for the exact density 
functional [Eq. (0)]) will be restored: 



ElDaIPq + Pj] + Ecorr[Pj] ~ ElDA 



+ / rfrp,(r)($(r) +\ j rfr'pr^(r')W^(r, r') 

J ^ Pj (r)=0 2 J 



Pj(r)=0 
Pi(r)=0> 



(10) 



Pj(r)=0 

Equation (fTUI) is equivalent to Eq. (0) with the function Aj(r) equal to 
A^(r) = $(r) _ ^ + ^ / rfr'pf"''(r')W^(r, r' 



p,(r)=0 2 . 

We define the "orbital densities functional" Eqdf as: 



Pj(r)=0 



(11) 



= i5;z.z.^[p] - ^ E / ^^/'^(r) / - Pr^(r ))W^(r, r') , • (12) 

^^J J ft(r)=0 

Please note that the ODE functional [Eq. (jl2|l ] depends not only on the total charge density 
p(r) that is a sum of orbital densities but also explicitly on a set of orbital densities {pj(r)}. 

Let us investigate the dependence of the functional [Eq. (fT^] on the number of electrons 
A^. At first please note that for the integer values of A^ corresponding to the integer orbital 
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occupancies rij the value of the functional [Eq. p2j)] coincides with the corresponding value 
of the LDA functional. For integer Ui values orbital density Pi{r) is equal either zero for 
empty state (n^ = 0) or p^°'^{r) for occupied state (n^ = 1). In both cases correction term 
in Eq. (jl2|) vanishes. We will show that in accordance with the properties^i of the the exact 
density functional Eedf [Eq- (0)] this dependence corresponds to the curve of Eq d p versus 

as a series of straight-line segments with derivative discontinuities at integer values of A^. 

The increase of the total number of electrons occurs via consequential increase of 
the orbital occupancies rij so that when A^ changes from M to M + 1 the value of the 
corresponding to the highest occupied orbital for the system with M + 1 electrons changes 
from to 1. The corresponding variation 5p(r) of the total charge density p(r) will consist 
exclusively of the variation of corresponding orbital density Pj(r) of particular orbital j. 
Then variational derivative of the functional [Eq. (fT^ ] is equal to: 

SEqof _ 5Eoop _ 6E1JXA [ ^ 1 ^maxl„l\ 

(5p(r) 5pj{v) (5p(r) 
The first term in the right part of Eq. (jl3p 6El£,a/Sp{t^) is continuous. However, the second 
term depends on the index j. When number of electrons A^ is equal to M — u then j 
corresponds to the highest occupied orbital for the system with M electrons. However, for 
N = M + uj the index will change to j ' + 1 corresponding to the lowest unoccupied orbital 
for the system with M electrons or the highest occupied orbital for the system with M + 1 
electrons. The value of Pj{r') in the integral in the right part of Eq. (jl3|) will jump from 
pmax^j.;^ to zero with A^ going M — uj to M + u for infinitesimally small u. That results in 
the corresponding jump of variational derivative 6Eodf/Sp{t^)'- 



j dv\p,{v')-\pJ^-{v'))W{v,v') . (13) 

J 2 Pi(r)=0 



SEqdf 



SEqdf 



5p{r) M+Lo 5p{r 
- I dY'(p'^''^{r')W{r,v') +p'rf(r')W^(r,r' 

2 J \ ■> Pj{r)=0 



M-IaJ 



Pj+i(r)=0/ 

In order to show that the variational derivative 5Eodf/ ^p{^) is constant for number of 
electrons varying from M to M + 1 we need explicit expression for the variational derivative 
of LDA functional 5Elda/ ^pij)- That could be obtained using expansion of Eqs. ((Tj) and 
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6p{r) 5pj{r) 5p{r) p,(r)=o J 5p{r)6p{r' 



i-'pM)^j§^. (15) 



p,(r)=0 



Mr] 



p,(r)=0 



(ir'pj(r')iy(r, r') 



Pj(r)=0 



Then Eq. fll3p takes the form: 

SEoDF SEoDF SElda 



Sp{r) 6p,{v) ' Mr) Mr)=o 2j ^ ' 6p{r)6p{r' 



(16) 

Pi (r)=0 



_^ + ^J dv'pf^^iv')Wir,r') 



Pj(r)=0 



The right part of Eq. ()16|) does not depend on pj and hence it is constant for fractional 
number of electrons rij. As derivatives of all terms in LDA functional are continuous the 
jump of SEodf/ ^p{^) in Eq- f!14p can be assigned to the jump in the variational derivative of 
the effective exchange-correlation energy term in ODF functional 5E'^^^ / 5p{v). This term 
can be derived from ODF functional [Eq. ()12j)] similarly to the LDA case Eq. ()A2|) . 

If orbital ipi{r) is fixed then variation of orbital density pj(r) occurs only via variation of 
occupancy rij: 

<5p,(r) = 5n#,(r)p. (17) 
Then the second variation of LDA functional is: 

5''Elda = \llj ^^^Z'^^^) / dv'5p,{v')W{v,v') (18) 



\Y.6nM, I dr\Uv)\'l dr'\i^,{r')\'W{r,r'). 
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The second derivative of LDA functional with respect to the occupancy rii is: 

^^^^ - I dr\Ur)? j dv'\i^,{v')\^W{vy). (19) 



5^Er 



5ni5n,j 

The value of the derivative 5'^Eioa/ 5ni5nj could be obtained in constrained LDA calcu- 
lation either in a direct way or by taking into account that derivative of LDA functional over 
orbital occupancy 5Elda/ is equal to one-electron eigenvalue [see Eqs. ()A26|) - ()A31|) ]: 

Using Eqs. |T7j) - pn|) ODF functional [Eq. |T^ ] could be rewritten via occupancies n^: 

EodfW), {n,}] = EldaW)] - ^ 5^n,(n, - 1)^. (21) 

i * 

The derivative of ODF functional [Eq. p2|) ] with respect to the total number of electrons 
can be expressed via occupancies Ui of partially occupied orbital j: 



SEqdf SEqdf SElda , /I ^Sej 1 6ej 

— ^ = — = \- — - = e,- + 

5N Snj 6nj ^2 ^'6nj ^ ^2 ^'Snj 



(22) 



Expanding e function of Uj and keeping only the first term (which is equivalent to 

keeping only first and second variational derivatives in Eq. ((7j)) one can show that derivative 
5Eodf/^N does not depend on nf. 



6E, 



ODF 

~ c 

' Co 



nj=0 orij 



(23) 



Analogue of Eq. (jl4p demonstrates a jump of bEoDpj when number of electrons N is 





1 btj 






nj =o 2 (5r;,j 


nj=0 



going through an integer value M: 



(5E 



+ + (24) 



^ "'"'""^ "'^ ' 2\5nj ' 6nj^i 



We have shown that with the accuracy of expansion of Eq. (|7j) the variational derivative 
of the "orbital densities functional" Eqdf [Eq. ()12|) ] conforms to the conditions for the exact 
density functional found in Ref. it is a constant for fractional number of electrons and 
has a discontinuity when number of electrons is going through an integer value. 

Equations and are directly related to "transition state" approach proposed 
by Slater— to calculate excitation energies. In this scheme LDA eigenvalue (Kohn-Sham 
equations eigenvalue) of the corresponding one-electron state should be calculated with its 
occupancy equal 0.5 (half way between initial and final states of excitation process). One can 
identify a derivative of the ODF functional [Eq. (j2H) ] over occupancy rij as a corresponding 
one-electron energy eP"^^: 

If LDA eigenvalue ej(nj) is a linear function of occupancy n.j {6ei/6ni = const) (which is 
equivalent to expansion Eq. ^I^i ) then for an empty state (n^ = 0): 

6f^^(0.5) = ef^^iO) + Ig. (26) 

For occupied states (n^ = 1) a sign plus in Eq. will be replaced by minus. One can see 
that Eq. can reproduce both those cases: 

e?^' = e^^\n.) + {\-n.)^^ = 6.^^^(0.5). (27) 



In Ref. |16| we have introduced an auxiliary functional by variation of which Eq. (j27|) giving 
"transition state" correction to one-electron energies can be obtained. That functional is 
identical to the functional [Eq. ()21|)]. Keeping in mind that we have defined functional 



[Eq. ()21|) ] in order to reproduce properties of the exact density functionalii one can arrive to 
conclusion that eigenvalues calculated in "transition state" method are a good approximation 
to the eigenvalues of the exact functional. 

The minimization of the functional Eq. (jl2|) will give a set of orbital densities {pj(r)} and 
hence the total charge density p(r) corresponding to the ground state of the system. It is 
useful to derive equations for fluctuations of the orbital densities {5pj(r)} around the average 
ground state functions {pj(r)}. This "ground state" charge density p(r) corresponds to the 
minimum of ODF functional [Eq. (fT^ ] but not of the LDA functional so that one must use 
in the calculations constrain potential $(r) (see Appendix |3 Eq. ()A25|) ). Analogous to the 
expansion Eqs. ((Tj) and (jH)) one can write for the LDA part of ODF functional: 

ELDA[{p^ + Sp,]] ^ ELDA[{Pi}] + [ dr6pi{r)^r) (28) 
+ ^EE/ ^^'^Z^^^^) / dr'6p,{r')W{r,T'). 

i j 

For the correction term Eq. the corresponding expression is: 

Ecorr [{p^ + S p,}] ^ E.^rr [{p^}] ' (29) 

dr6p.{r) I dr'(p,(r')-ipr^(r'))W^(r,r') 

Both Eqs. and become exact if one supposes that /p(^^)Sp{r') ~ ^(^' '^'^^^ 
depend on fluctuations Spi{r) so that all variational derivatives higher than the second order 
are equal to zero. In the following we assume that this approximation is valid. Then for the 
total ODF functional one has: 

EoDF[{pi + Spi}] = EoDF[{pi}] + (30) 

+ drSp,{r)(^^r)- J rfr'(p,(r') - ipr^(r'))W^(r, r')) 

+ / dv'5p,{v')Wir,v'). 

It is convenient to introduce effective ODF potential Vf^^{T): 

\/«^^(r) = $(r) - J dv\p,{v') - ^pr^(r'))IV(r, r'). (31) 
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Using Eqs. ()30p and ()31|) one can define a Hamiltonian for tfie density matrix fluctuation 
operators 5pj(r) = Pi{r) — Pi{r) (here ground state orbital density can be considered as an 
average value of density matrix operator Pi(r) = (pi(r))): 

i i j^i 

The orbital density matrix operators Pi(r) could be expressed via orbital occupancy op- 
erators rii (if orbitals il)i{Y) are fixed and the variation of orbital density pj(r) occurs only 
via variation of occupancies rii): 

U^)=n,\^,{v)\\ (33) 

The fluctuation Hamiltonian Eq. (jH^ can be expressed using Eqs. ((T7j)-(j2n)) via occupancy 
operators and their average values rij = (nj) [drii =ni — {%)) 

H = J2ie. + (1 - n.)^)6;^. + ^ E E (^4) 

i ^ i j^i ^ 

The orbital densities Pi(r) are defined by orbitals ipi{r) in expression for charge den- 
sity in Eq. In the density functional theory it is usually assumed that those orbitals 
are solutions of Kohn-Sham equations [Eqs. ()A3|) and ()A15|) ] and so minimize functional 
Eq. ()A2j) . However, any unitary transformation (defined by unitary matrix U) of the set of 
the functions ipi{r) produces a new set of orbitals: 

(35) 

j 

corresponding to the same charge density p(r). This new set of orbitals can be used to define 
ODF functional in Eq. p2|) an so this functional is orbital-dependent. 

In order to remove uncertainty in choosing the orbital set one needs to impose an addi- 
tional condition. We propose the following way to do it. The Hamiltonian Eq. (jH^ can be 
used to calculate the contribution to the energy from fluctuations: 

= ^ E E / / dr'{6~p,{r)6~p^{v'))Wir, r'). (36) 

(The average value of fluctuations (5pj(r)) is equal to zero so the first term in Eq. (jS^I) does 
not give contribution to the fluctuation energy.) Taking into account Eqs. (f^Hj) and (f^^ one 
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can separate contributions to Eq. ()36|) from LDA functional where summation over i and j 
is performed including terms with i = j and a correction term Eq. (j2H). 
The LDA part is: 

(Hlda) = lj2J2jd^J dr'{6^pAr)6^p^{r'))W{r,r') (37) 

i j 

= ljdrj rfr'(5^4(r)5^5;o,.(r'))l^(r,r') 



dr j dr'{5p{r)5p{r'))W{v,r'). 
Correction term gives a negative contribution equal to: 



{H^orr) = E / / dr'{S'p,{r)fp,{r'))W{r, r'). (38) 

From Eq. ()37p one can see that LDA contribution to fluctuation energy is defined by 
the total charge density fluctuations {6p{r)6p{r')) and so does not depend on the orbitals 
definition. Then the minimum of the fluctuation energy Eq. ()3fj|l is achieved when the 
absolute value of correction contribution Eq. has a maximum. 

The correction term contribution in Eq. ()38j) can be calculated using expression of fluc- 
tuations Hamiltonian via fluctuation occupancy operators 6ni [Eq. (jH^ J: 

^ ]^ . Q^. 

{Hcorr) = -- (39) 
i * 

As the average value of square of occupancy fluctuations {6ni6ni) depends on the specific 
properties of the system, the only way to minimize fluctuation energy (H) [Eq. ()36j) ] is to 
maximize dei/dui. Using Eqs. and (pn|l this parameter can be expressed via orbitals 
as: 



de, S^E 



For new set of the functions V'i(r) [Eq. (j35|l ] one has 



''''^ - j dv\i;,{v)\^ j dv'\i;,{v')\^W{vy). (40) 



^ = ^ f/,,t/*.,f/,,f/*, I dv I rfrV,(r)^;,(r)l^(r,r')V^Kr')^r'(r')- (41) 

Using Eq. ()4H1 and the fact that the derivatives beijbni are always positive one can define 
a functional of unitary matrix \J : 

^[^] = E = E E ^^^l^^^l^ [ c/rV,(r)^;,(r)W^(r, r')^Kr')#(r') (42) 



i jj'll 
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maximization of which one can use as a condition to determine the matrix U and hence the 
optimal set of orbitals ipi{r) to define ODF functionaL 

The function W{r, r') [Eq. ()A24|) ] is defined as screened effective interaction between 
density fluctuations Sp{r), Sp{r') and so should decay with increasing of |r — r'| value. Then 
the more localized in space orbitals ipi{r) are, the larger should be the integral Eq. (|4(jp 
value. One of possible choices for these orbitals could be Wannier functions. Marzari 
and Vanderbilt in Ref. Q proposed the condition of maximum localization to determine 
the procedure to calculate Wannier functions for multi-band case (two orthonormal sets 
of functions, Wannier and Bloch, are connected via unitary transformation so Wannier 
functions can be considered as a particular choice of unitary matrix U in Eq. (|35p). The 
requirement of maximum localization should lead to the reasonably maximized values of 
parameter dei/drii in Eq. (pn|) and hence to minimization the fluctuations energy Eq. (jHUjl . 
Then Wannier functions obtained via procedure proposed in Ref. |l5| are a good choice for a 
set of orbitals to define ODF functional. 

In the following we will assume that orbital densities are defined not by one-electron 
eigenfunctions ipi{r) as in Eq. (jH) but by Wannier functions Wn{r) and their occupancies Qn 
calculated via procedure described in Sec. IIIII 

p(r) = 5^p„(r), (43) 

n 

Pn{r)=Qn\Wn{r)\'. 

Using Eqs. (PT|) -()33 p ODF functional [Eq. (fT^ J could be rewritten with Wannier functions 
occupancy operators Qn [Eq. ()B1|) ] and their average values Qn = {Qn)- The corresponding 
functional will be analogous to Eq. (PT|) but with occupancies Qn [Eq. and energies En 
[Eq. (jnni)] corresponding to Wannier functions: 

1 dE 

EoDF = Elba - - j^QniQn - !)• (44) 

The variation of the functional Eq. ()44|) will produce (see Appendix El one-electron 
Hamiltonian Hq^p in the form of projection operator: 

H'ODF = HlDA + = ^LD^ + H \^nWn{W^l (45) 

n nT 

\W'^) in Eq. are Wannier functions [Eq. and 5Vn are: 
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The values of derivatives dEn/dQn should be determined in constrained LDA calculations 
(see Seclini). 

Equations are identical to the equations of "generalized transition state" (GTS) 

method that we have developed in Ref. Q. The basis of GTS method was an idea that one- 
electron energies corresponding to Wannier functions should have a meaning of the removal 
(addition) energies for electrons from (to) the corresponding states. That was realized by 
using "transition state" scheme^^ generalized to Wannier functions instead of eigenfunctions. 
It is remarkable that the same equations can be obtained by introducing correction to the 
LDA functional restoring the properties of the exact density functional theory.— 

The potential correction operator [Eqs. (|45|) and (j46p ] shifts the energies of Wannier func- 
tions on 5Vn values which are negative {5Vn = ^||^) for the occupied states {Qn = 1) 
and positive {5Vn = |f^) for the empty states {Qn = 0). That gives valence bands pushed 
down and conduction bands pushed up comparing with standard LDA results increasing en- 
ergy separation between them (energy gap value) which is systematically underestimated in 
LDA. It was shown in Ref. Q that calculations with potential correction [Eqs. ()45|1 and (jlEI)] 
result in a much better agreement with experimental energy gap values for semiconductor 
(Si), band insulator (MgO), Mott insulator (NiO), and Peierls insulator BaBiOs. 

The Hamiltonian H^^p [Eq. ^] can be considered as a static mean-field approximation 
to a general problem. Wannier functions analogue of the fiuctuation Hamiltonian Eq. (jM)) 
is: 

HODF = H'oDF + 2 E E - {Qn)){Qn' " (Qn'))- (47) 

The first and second parts of Hamiltonian Eq. (j47|) are not (!) noninteracting Hamiltonian 
and interaction term as it is usually defined in model Hubbard and Anderson Hamiltoni- 
ans. The first part -f^oDF is equivalent to the Hartree-Fock approximation Hamiltonian 
determined by the average values of Wannier functions occupancies {Qn) and the second 
part describes interaction between fiuctuations around (Qn)- As these average values are 
determined from the solution of the full Hamiltonian Eq. (|47|) . that defines a self-consistent 
calculation scheme. In contrast to LDA+U^ and LDA+DMFT— methods there is no 
"double counting" problem in this Hamiltonian because there were no "merging" of LDA 
and Hubbard model concepts here and both terms in Eq. ^7\\ were derived from the same 
functional Eq. dHj). 
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The problem defined by the Hamiltonian Eq. ()47p can be solved by any of the methods 
developed to treat many-body effects. In the present work we have used dynamical mean- 
field theory (DMFT)i2JL3*i4 ^^^^ ggp_ ^ 

III. DEFINITION AND CONSTRUCTION OF WANNIER FUNCTIONS 

The orbital projection calculation scheme for Wannier functions (WFs) used in the present 
work was described in details in the earlier papev^ where the LDA+DMFT (DMFT - 
dynamical mean-field theory) method in the Wannier function basis set was proposed. Below 
we present the main formulas of this scheme. 

The concept of WFs has a very important place in the electron theory in solids since its 
first introduction in 1937 by Wannier.— WFs are the Fourier transformation of Bloch states 

IV'ik), 

\W^) = ^^E^"^'"!^^-^)' (48) 



where is the number of discrete k points in the first Brillouin zone (or, the number of 
cells in the crystal) and T is lattice translation vector. 

Wannier functions are not uniquely defined for a many-band case because for a certain set 
of bands any orthogonal linear combination of Bloch functions \ipi-k) can be used in Eq. 



In general it means the freedom to choose the of unitary transformation matrix U^^^ for 
Bloch functions)^ 

l^^k) = Ef^i'V.k). (49) 

j 

There is no rigorous way to define Ujf\ This calls for an additional restriction on the 
properties of WFs. Among others, Marzari and Vanderbillii^ used the condition of maximum 
localization for WFs, resulting in a variational procedure to calculate Ujf\ To get a good 
initial guess the authors of Ref. Q proposed choosing a set of localized trial orbitals 
and projecting them onto the Bloch functions It was found that this starting guess is 

usually quite good. This fact later led to the simplified calculating scheme in Ref. |2l| where 
the variational procedure was abandoned and the result of the aforementioned projection 
was considered as the final step. 
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For the projection procedure used in the present work one needs to identify the set of 
bands and corresponding set of locahzed trial orbitals |0n)- The choice of bands and orbitals 
is determined by the physics of the system under consideration. 

The set of bands can be defined either by the band indices of the corresponding Bloch 
functions {Ni, N2), or by choosing the energy interval {Ei,E2) in which the bands are 
located. Nonorthogonalized WFs in reciprocal space |VF,ik) are then the projection of the 
set of site-centered atomiclike trial orbitals on the Bloch functions of th chosen 
bands defined by band indices (A^^i to N2) or by energy interval {Ei,E2): 

\Wnl,) = = Yl l^^k)(^^k|0n). (50) 

i=Ni i{Ei<ei('k)<E2) 

In the present work we have used LMTO method^^ to solve a band structure problem 
and the trial orbitals were LMTOs. The coefficients (^/'jkl^n) in Eq- fj^Hj) define (after 
orthonormalization) the unitary transformation matrix f/jf^ in Eq. ()49|1 . 

In any DFT method the Kohn-Sham orbitals are expanded through the certain basis: 

The basis functions of the LMTO method are Bloch sums of the site-centered orbitals: 



J=5^e''^^0,(r-R,-T), (52) 



where /i is the combined index representing qlm [q is the atomic number index in the unit 
cell, Im are orbital and magnetic quantum numbers), Rg is the position of atom in the unit 
cell. 

For the orthogonal LMTO basis cj^j = {(f>ij.\ipik.) and hence 



^ N2 N2 

i=Ni i=Ni fi fj. 



with 

N2 

bin = E (54) 

i=Ni 

For a nonorthogonal basis set orthogonalization of the Hamiltonian must be done before 
using Eq. (|^. 
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In order to orthonormalize the WFs [Eq. ()53|) ] one needs to calculate the overlap matrix 
On„'(k), 

and its inverse square root 5'„„/(k), 

Snn\\^) = Oj'(k). (56) 

From Eqs. ()53p and ()56p the orthonormalized WFs in k space |VF„k) can be obtained as 
iW^nk) = $^5„.,(k)|W^„,k) = I^^i')^™ = E^nl'^J)' (57) 

n' i=A''i /J 



with 



4* = (^.kllVnk) = E ^nn'(k)c!:r„ (58) 



= (0j|W^nk) = 5^ Cj,4*- (59) 



Using Eqs. one can find energies of WFs: 



= j^T.T.'^n&m (60) 



k i=Afi 

and their occupancies: 

Qn = (W^„^l(X;iV^.k)^(i?/-e.(k))(^,k|)|W^„^) 

N2 

= J;^Y.Y.^n^Cn'AEf-em, (6I) 

k i=Ni 

where ej(k) is the eigenvalue for a particular band, 6{x) is the step function, Ej is the Fermi 
energy. 

The transformation from LMTO to the WF basis set is defined by the explicit form of 
WFs Eqs. ()57|) and (|^^. and by the expressions for matrix elements of the Hamiltonian 
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and density matrix operators in WF basis [Eqs. ()6Up and (jHU)]- The transformation from 
WF to LMTO basis can also be defined using Eq. (jFfj) . Such transformation is necessary 
in calculations using ODF correction potential in the form of projection operator Eq. (|45|) . 
It will also be used for constrained LDA calculations needed to determine the values of 
WFs energy derivatives with respect to their occupancies dEn/dQn' which define potential 
correction parameters 5Vn in Eqs. fj45|) and PU)) and effective Coulomb parameters Unn' in 
ODF fluctuation Hamiltonian in the Hubbard form Eqs. (f7Bj) - (f7Sj) . For example if constrain 
potential operator is diagonal in WF basis = A„(5„„/), then in LMTO basis its matrix 

elements can be calculated via: 

^const.^^|^T^^^^^T|^ (62) 



\W^) = Y.^-'''''h)^\^)). (63) 
//--*^(k) = = E KnKbll (64) 



IV. WANNIER FUNCTIONS IN THE GREEN-FUNCTION FORMALISM 

In many-body theory the system is usually described not by Bloch functions |'?/'ik) 
[Eq. (j^H)] and their energies ej(k) but by the Green function 

^(r' = ^ E ^'(r' ^) = ^ E E <^S(r)Gj.(5)0:^r'). (65) 

k k i-iv 

The matrix Green function G^j^{e) is defined via the noninter acting Hamiltonian matrix 
if°^(k) [Eq. dZHl) and ^] and the matrix self-energy S}^^(5) [Eq. ^] as 

Gl,{e) = [e- H\k) - k) + zr/];J. (66) 

We define the nonorthonormalized WF obtained by projecting the trial orbital (/'n(r) on the 
Hilbert subspace determined by the Green function [Eq. ()65p ] in the energy interval {Ei, E2), 
namely, 

WnUr) = — Im / de dr'G\v, r', e)^^{v') = ^ &J„0j(r) (67) 
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and 



1 



TT 



E2 



^ --Im / deGUe). (68) 



El 



In the noninteracting case, the self-energy S(e, k) is absent and hence we have 

where cj^j are the eigenvectors, and ei(k) are the eigenvalues of H^{k). Thus in Eq. (|HH)l 
becomes 

N2 

i=Ni 

where Ni, N2 are the band numbers corresponding to the energy interval {Ei, E2). Since this 
recovers the result of Eq. (jSljl . we demonstrated that our general definition of WFs [Eq. (|67j) ] 
via Green function reduces to that in terms of Bloch functions [Eq. (I53|) ] in Sec. Illll 

To orthonormalize IVnk(r) defined in Eq. (j67|) . one can just follow the orthonormalizing 
procedure made in Sec. IIIII [Eqs. (j^ -(fS ^ ] with the overlap matrix 0„„/(k) is defined as 

The occupancy matrix in the orthogonalized WF basis [Eq. (|U7|)] is defined as 

Q„„,(T) = -hm p de j j drdr'^ ^ iy4(r)G'''(r, r', £)iy„,k(r')e-'''^. (71) 
By using Eq. (jHS)) and orthogonalized Wannier functions [Eq. fl67j) ]. one finds 

Qnn'(T) = ^ E E ^>ln'Q>~'''^ . (72) 

with 

Ql, = --Im r'dsGl^ie). (73) 

^ J — oo 



The energy matrix can be defined similarly (except that the integral over energy is cal- 
culated in the {Ei,E2) interval where the corresponding WFs are determined) as 

-Ei 



--Im rede f I drdr'l ^ ^nk(r)G''(r, r', £)iy„,k(r')e-^''^ (74) 



1 \^\^Lk*7k pk ^-ikT 

k iiu 



19 



with 



1. 



TT 



E2 



Ej, = — Im / edeGUe). (75) 



El 



V. DMFT FOR ODF FLUCTUATION HAMILTONIAN 

The DMFTii^ii^iM was recently found to be a powerful tool to numerically solve multiband 
Hubbard models. In order to use this tool the fluctuation Hamiltonian [Eq. ()47|) ] should 
be rewritten in the form of standard multi-orbital Hubbard model. For that one needs to 
identify Coulomb parameters f/„„/ as derivatives dEn/dQn' and rearrange terms in Eq. (j47p 
into noninteracting and interaction parts: 

HoDF = + Hint (76) 

with interaction term: 

Hint = ^ Unn'QnQn' (77) 

and noninteracting Hamiltonian: 

H^ = Hq^p + ^ Unn'{Qn){Qn') ~ Qn Unn'jQn')- (78) 

In DMFT the lattice problem becomes an effective single-site problem which has to be 
solved self-consistently for the matrix self energy S and the local matrix Green function in 
Wannier functions basis set: 

Gnn'ie)=-^ [ dk(\ie-fx)l-H'ik)-^ie)Y'] , (79) 

VbZ J V L ^ / nn' 

where yU is a chemical potential, H^ is noninteracting Hamiltonian [Eq. ()78|)] and Ti{e) is 
self-energy in Wannier functions basis: 

E{e) = J2\Wn)^nn'{e){Wn'\. (80) 

nn' 

The DMFT single-site problem may be viewed as a self-consistent single-impurity Ander- 
son model.— The corresponding local one-particle matrix Green function G can be written as 
a functional integraU^ involving an action where the Hamiltonian of the correlation problem 



20 



under investigation, including the interaction term with the Hubbard interaction, enters.— 
The action depends on the bath matrix Green function Q through 



-1 



(gf =(G)"i + S. (81) 

To solve the functional integral of the effective single-impurity Anderson problem, various 
methods can be used: quantum Monte Carlo (QMC), numerical renormalization group 
(NRG), exact diagonalization (ED), non-crossing approximation (NCA), etc. (for a brief 



overview of the methods see Ref. 
solve the impurity problem. 



18^ . In the present work the QMC method was used to 



VI. SELF-CONSISTENCY 

The matrix self-energy Tj{e) obtained as a solution of DMFT in Sec. |V]is defined in the 
WF basis set [Eq. fl80|) ]. In order to compute the interacting Green function in the full- 
orbital Hilbert space [Eqs. ()65|) and ()66|) ] one has to transform self-energy to the LMTO 
basis set. This can be easily done by using the linear-expansion form of the WFs in terms 
of the LMTO basis set [Eqs. (jHTj) and (I3UD ]. 

E;^,(£) = (0j|S(5)|0^) = 5^(0j|iy.k)Snn'(^)(W^n'k|0^) = J] (e)^, . (82) 

nn' nn' 

The matrix elements of the self-energy Sj^^(£:) [Eq. ()82j)] together with the noninteracting 
Hamiltonian matrix iJ|]^(k) [Eqs. ()78|) and (jMj) ] allow one to calculate the matrix Green 
function G^^,{e) [Eq. and thus the interacting Green function G{r,r',e) [Eq. 

G{r, r', e) contains the full information about the system, and various electronic, magnetic, 
and spectral properties can be obtained from it. 

One can also calculate the charge-density distribution modified by correlation effects via 

p(r) = --Im / ^ deG{r,r,e). (83) 

With this p(r) one can recalculate the LDA potential which is a functional of electron density 
[Eq. ()A7|) ]. From the Green function [Eq. (j65|) ] one can recalculate new WFs [Eqs. (j67|) and 
(|68p] which together with the new LDA Hamiltonian allows one to obtain new parameters 
for the noninteracting Hamiltonian [Eqs. (jMj) and ()78p]. With the new set of Wannier func- 
tions one performs a series of constrain LDA calculation to determine derivatives dEn/dQn' 

21 



and hence define a new Coulomb interaction parameters Unn' in the interaction Hamilto- 
nian [Eq. (f77j) ]. The set of new LDA potential, WFs, and Coulomb interaction parameters 
calculated from the interacting Green function [Eq. (j65|) ] defines the input for the next it- 
eration step and hence closes the self-consistency loop in the computation scheme. For the 
feedback from DMFT to LDA in the particular case of the LMTO method^^ one needs a set 
of moments for the partial densities of states Alj^f for every atomic sphere q and the orbital 
moment I (Ref. |23|) in order to calculate the new charge density and hence the new LDA 



potential: 



M^f = I ' deePNgiie), (84) 



k m 

The orbital densities functional computational scheme described above is ab initio (does 
not contain any outside parameters) and fully self-consistent. The charge density p(r), Wan- 
nier functions Wn{r), and derivatives dEn/dQn' (effective fluctuations interaction strength) 
are recalculated on every self-consistency loop and hence modified by correlations comparing 
with the values obtained in standard LDA calculations. Analogous calculation scheme was 
proposed in Fig. 1 in Ref. Q. 

In the present work ODF method was applied to the problem of metal-insulator tran- 
sition in nonstoichiometric lanthanum trihydride LaH3„^., where both parts of the calcula- 
tion scheme, static mean-field approximation [Eqs. ()44|l - PH|l ] and fluctuation Hamiltonian 
[Eq. ()47|l ] in DMFT (Sec. |^ were needed to describe experimentally observed dependence 
of the system ground state on the hydrogen deficiency parameter x. 



VII. ELECTRONIC STRUCTURE OF NONSTOICHIOMETRIC LANTHANUM 
TRIHYDRIDE 

Lanthanum trihydride LaH3_^ shows an interesting metal-insulator transition with in- 
creasing of X value.— While stoichiometric LaHa is an insulator, 30% hydrogen deficiency 
(x ~ 0.3) results in a metallic ground state. However, standard LDA calculations^^ give a 
metal as a ground state even for stoichiometric composition LaHs. In this case well-known 
problem of LDA, underestimation of the energy gap value, is so severe that the gap value is 
negative with valence and conduction bands overlapping on 0.25 eV (see Fig. [TJ. 
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The lanthanum trihydride has a crystal structure derived from the face-centered cubic 
structure (lattice parameter a = 10.5946 a.u.) for La atoms and hydrogen atoms occupying 
two tetrahedral and one octahedral interstitials per every metallic ion. The band structure 
of LaHs (Fig. Q]) is relatively simple: lower three occupied bands are formed by hydrogen Is 
states and the conduction bands correspond to lanthanum 5d, 6s, and 6p states. 

Attempts to cure the LDA fault in LaHa were performed with many methods beyond 
LDA, among others, the GW method^^ (the energy gap Eg=0.8-0.9 eV), the weighted local 
density approach^! (^Eg=0.7 eV), and model calculations.— 

The orbital densities functional (ODF) theory [Eq. ()44|l ] proposed in the present work 
adds to standard Kohn-Sham equations potential correction in the form of projection opera- 
tor [Eqs. (j45p and (j46|) ]. This correction is negative for occupied valence bands and positive 
for empty conduction bands and hence increases energy separation between those bands. 
ODF calculations for LaHs results in insulating ground state with a fundamental gap value 
of 1.10 eV and direct optical gap 1.25 eV at G-point (see Fig. |2)). Experimental data for 
the value ot energy gap E,^0.5 eV wa. e.«„>ated u.ng act.vat.on energy deternnned fron. 
resistivity measurements in Ref. |22l Optical measurements estimate the direct gap E^ as 
1.87 eV but the fundamental band gap ~1 eV lower.— 

Hydrogen atoms could be removed from the lanthanum trihydride forming nonstoichio- 
metric LaH3_a;. We have performed LDA and ODF calculations for two compositions x=0.25 
(LaH2.75) (see Fig. 01) and x=0.5 (LaHg.s) (see Fig.|H). 

Both LDA and ODF calculations gave metallic ground state for two compositions x=0.25 
and x=0.5. Every hydrogen vacancy leads to appearance of one vacancy state in the gap 
split from the conduction band with one electron occupying this state (see Fig. Ofa)). Static 
mean-field ODF potential correction [Eqs. (^3)) and ()46p ] results in separation of the vacancy 
band from the conduction band and in increasing of energy gap between valence and conduc- 
tion bands comparing with LDA but half-filled vacancy band stays metallic (see Figs. El^b) 
andllfb)). However, experimentally only LaH2.5 is metallic while LaH2.75 should be an insu- 
lator. To treat this problem we have used ODF fluctuation Hamiltonian [Eq. (@Z|)] and have 
solved it using DMFT-QMC (see Sec. |V}. Inverse temperature parameter was /5=10 eV~^. 
In the result we have obtained in agreement with experiment paramagnetic insulator ground 
state for LaH2.75 with a typical Mott insulator pattern of lower and upper Hubbard bands 
around chemical potential with the energy gap of 0.1 eV (see Fig.|31^b)). However, for LaH2.5 
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a width of the vacancy band is significantly larger then for LaH2.75 (see Figs. El^b) and|3fb)) 
and DMFT calculations resulted in well defined metallic state (see Fig. Efa)). 



VIII. CONCLUSION 

We have proposed orbital densities functional that by construction has discontinuous 
exchange-correlation potential in agreement with exact density functional theory. The one- 
electron potential obtained by variation of this functional produces comparing with LDA 
lower energies for valence bands and higher energies for conduction bands thus overcoming 
systematic underestimation of energy gap value in LDA. To treat correlation effects we have 
derived from orbital densities functional the Hamiltonian corresponding to the fluctuations 
of orbital densities around the ground state values. Combining this Hamiltonian with dy- 
namical mean-flled theory and quantum Monte Carlo method for effective impurity problem 
we have developed ab initio and fully self-consistent scheme for electronic structure calcula- 
tions. The scheme was applied to the problem of metal-insulator transition in LaHs-^: and 
gave a qualitative improvement for agreement with experimental data. 
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APPENDIX A: CONSTRAIN LDA FUNCTIONAL AND ITS DERIVATIVES 

LDA functional Eloa[p] depends on the electron density p(r) deflned by one-electron 
functions ipi{r) and their occupancies nf. 




(Al) 
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Elda[p] = T[p] + ^Jdr Jdr' 




(A2) 



+ / drVz{r)p{r) + E,,[p], 



where Vz(r) is external potential, Exc[p] is exchange-correlation energy functional [see 
Eq. (^] and T[p] is kinetic energy term. 

Kinetic energy is a functional of electron density T[p] via single particle orbitals ipi{r) 
which are functionals of density themselves. To prove that the following reasoning can be 
used.— Hohenberg-Kohn theoren*^ states that external potential is uniquely defined by the 
ground state electron density of the system. That is also true for a noninteracting system 
with potential V5[p](r). This potential defines a Schrodinger equation for the orbitals ipi{r): 



As potential V5[p](r) is a functional of density then solutions of the Eq. ()A3|) are also 
functionals of density: 



Kinetic energy is defined by a set of orbitals ipi[p]{T) and hence it is also a functional of 
density: 



(A3) 



ipi{r) = ipi[p]{r). 



(A4) 




(A5) 



Minimization of the functional Eq. ()A2|) gives: 

6Elda[p] _ 6T[p] , ^ 



+ VLDA{r) = 



(A6) 



Sp{r) 6p{r) 
where one-electron LDA potential is: 




(AT) 



(Vxcir) = SExc[p]/Sp{r) is exchange-correlation potential). 
Variation of kinetic energy functional is: 




(A8) 
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First order perturbation theory for one-electron energies ei[p\ gives: 

5eAp] = J dr^|J:[p]{r){6Vs[p]{r))Mp]{r) (A9) 

Then Eq. (HH) is: 

5T[p] = -J dv'Vs[p]{r)5p{T) (AlO) 
and variational derivative of kinetic energy functional is: 

The minimization condition [Eq. ()A6|) ] becomes: 

^-^^ = -Vs[p]ir) + V,nA{r) = 0. (A12) 

A set of Kohn-Shan>i equations ()A3|) with y5[p](r) = VLDAi^) defines one-electron orbitals 
ipi[p]{r) and hence a new electron density p(r) via Eq. ()A1|) . The problem should be solved 
self-consistently till electron density used to calculate Vlda(i") will coincide with the density 
obtained using Kohn-Sham equations solutions ipi[p]{r). 

In order to define functional [Eq. ()A2|) ] for charge density p(r) different from the density 
minimizing Elda[p] one needs to introduce constrain potential $(r) and minimize auxiliary 
functional: 

F[p] = EldaIp] + J rfr$(r)(p(r) - p(r)), (A13) 

^ = + m = -Vs[p]{r) + V^nAir) + $ (r) = 0. (A14) 

Then potential in Kohn-Sham equations ()A3|) is: 



Vs[p]{r) = VLDA{r) + ^r). (A15) 

The second variational derivative of LDA functional is: 

6^Elda[p] S SVLDAjr) SVs[p]{r) 

mw) = M^^-'^f'^^'^ + "^""^'^^ = -jk^) — WV ^ 

The variational derivative of LDA potential [Eq. ()A7|) ] is obtained in a straightforward way: 

^^^^^(^) _ _J_ + - r')^ fAlT) 
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The second term in Eq. ()A16|) . SVs[p](r) / S p{r') is more complicated. Let us first solve an 
inverse problem - to calculate response function:— 

First-order perturbation theory for eigenf unctions and eigenvalues of Kohn-Sham equations 
(jX3|) gives: 

J^^=rdr'mr'). (A20) 
Then response function in Eq. ()A18|) is: 

The second term in Eq. ()A16|) is defined by the inverse response function: 
This inverse response function can be obtained solving integral equation: 



I dr"x{v,r")x^\r",r') = 6{r-r'). 



(A23) 



Then the second variational derivative of LDA functional is: 

^"^^LDA ^ 1 ^ _ ^/^ ^Kc(r) _ ^ ^^^^ ^^24) 

5p{r)6p{r') |r — r'| 6p{r') 

W{r, r') can be interpreted as an effective interaction strength between density fluctuations 

6p{r),6p{r') in LDA functional. The first variational derivative of the LDA functional can 

be obtained from Eq. ()A14jl : 

where $[p](r) is constrain potential needed to obtain electron density p(r) in self-consistent 
solution of Kohn-Sham equations [Eqs. ()A3j) and ()A15|l ]. 

Let us consider variation of the LDA functional Eq. ()A2|) in respect to the occupancies 
rij and orbitals ipi{r) defining electron density p(r) in Eq. ()A1|) : 

*^--E(*".^ + /«*wt§^). (A2«) 
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Srii Srii J ^p(r) 

rfrV;;(r)(-V^)^,(r)+ [ drrdr)VLDA{r)Ur) 



= n,{-V'mr) + n,V,nA{rmr) (A28) 
= mi-V' + VLDAir))^*{r) = n,e,^;(r), 

SElda = {^^^^i + ^^^i / c?r#i(r)V^*(r)). (A29) 
As orbitals are normalized: 

dr6ij,{r)tlj*{r) = 6{ [ rfr^,(r)^*(r) - 1) = (A30) 



and from Eq. (|A29|) it follows that: 

SElda 



e^. (A31) 



5ni 

All this was proven to be valid for the LDA eigenf unctions which are solutions of 

Kohn-Sham equations [Eq. ()A3|) ]. Let us prove that the derivative of LDA functional over 
the Wannier function occupancy Qn [Eq. (jH!))] is equal to the energy En [Eq. (p]j) ] of this 
WF. 

The charge density expressed via Wannier functions PVn(r) [Eq. ()57|) ] and its occupancies 
Qn is: 

/^w = EE^"i^"(r-T)r- (A32) 

T n 

Then derivative over Wannier function occupancy Qn (with fixed functions iy„(r)) can be 
calculated using Eq. ()A27|) : 



5Qn 



J drW:{T){-V^ + V{T))Wn{r) = J drW:{r)HLDAWn{r) = En- (A33) 
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APPENDIX B: PROJECTION OPERATOR IN ODF FUNCTIONAL 



In ODF functional Eq. pij) Wannier functions occupancies Qn [Eq. (jHH)] can be expressed 
via projection operator 

Qn = J2^HQn\i^^), (Bl) 

i 

(summation is over occupied states only in all equations below). The variational deriva- 
tive of occupancy Qn with respect to the one-electron wave function 6ipi{r) is: 

(B2) 

Then variation of the correction term in ODF functional Eq. ()44j) will give the following 
correction to the corresponding expression for LDA functional [see Eq. flA28|) ]: 

dEpDF _ ^^LDA_-^dEn I. 6Qn 

= (-V^ + KLDA(r))^;(r) + - Qn)Qn^PKr) 

= {-"S/"^ + VLDA{r) + J2^^nQn)^*{r). 

n 

This proves Eqs. (j45p and (j46p for ODF correction to LDA Hamiltonian. 



APPENDIX C: TOTAL ENERGY IN LDA AND ODF FUNCTIONALS 

Kinetic energy term in Eq. ()A2|) can be calculated using Eq. ()A3|) : 

(-V2)^i(r) = (e, - VLDAir))U^). (CI) 

Multiplying Eq. ()C1|) by ip*{r), integrating over r and summing over occupied states i one 
can obtain: 

J2 J dr^|J:{r){-V')Mr) = - /" rfr\4z,^(r)p(r). (C2) 

i i 

Then total energy in Eq. ()A2|) is equal to: 

Elba = Y.'^- j dvVLDA{r)p{r) + \ j j d^'^yz^ (C3) 
+ J drVz{r)p{r) + j dre^c{p{r))p{r). 
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Using Eq. ()A7|) for LDA potential V{r) the final expression for total energy is: 

Elda = " ^ / ^r/ ^"'iT^ + / ^r(^.e(p(r)) - \4e(p(r)))p(r). (C4) 

i 

In the case of ODF functional [Eq. (j44p ] there will be an additional term in expression 
for kinetic energy [Eq. ()C2|) ] due to the ODF potential correction [Eq. ()B3|) ]: 

drij*{r){Y,SynQn)Mr) = -J^^KQn. (C5) 

in n 

The additional terms to total energy expression Eq. ()C4|) are 

EoDF = Elda - J2 " 2 5Z 5Q^^"(^" " ^) (^^^ 

n n 

Here Elda is defined as right part of Eq. ()U4|) . -Elda differs from -E^da because it is taken 
for p(r) which does not minimize Elda- 
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FIG. 1: LaHs band structure and density of states calculated in standard LDA. 
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FIG. 2: LaHs band structure and density of states calculated in orbital densities functional (ODF) 
theory [Eq. (jUJ]. Zero of energy is at the Fermi energy. 
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FIG. 3: LaH2.75 density of states calculated (a) in standard LDA and (b) in orbital densities 
functional (ODF) theory [Eq. (gH)]. 
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FIG. 4: LaH2.5 density of states calculated (a) in standard LDA and (b) in orbital densities 
functional (ODF) theory [Eq. (gH)]. 
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FIG. 5: (a) LaH2.5 and (b) LaH2,75 densities of states calculated with ODF fluctuation Hamiltonian 
[Eq. (EZl)] in DMFT-QMC. 
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